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A  Functionally  Expanded  Computer  Program 
For  Computing  Machine-Dependent  Constants 


Introduction 

In  a  previous  paper  [1],  a  computer  program  was  discussed  that  can  be  used 
(independently  of  any  computer  vendor  claims)  as  a  benchmark  to  reveal  floating¬ 
point  arithmetic  characteristics  of  a  binary  device  computer  having  a  power  of  2 
base,  B(=  2k,  k  a  positive  integer),  and  in  transporting  mathematical  computer 
software  that  depend  on  the  machine’s  relative  accuracy  (computer  unit  roundoff 
error)  and/or  arithmetic  range  from  one  computer  to  another.  Furthermore,  since 
the  relative  accuracy  depends  on  the  minimum  (effective)  number  of  significant  bits 
in  the  mantissa  (fraction)  of  a  floating-point  number  and  the  roundoff  property  of 
the  computer  arithmetic,  the  program  also  computed  the  effective  number  of  bits  in 
the  mantissa  and  whether  the  arithmetic  rounds  or  chops.  However,  the  computer 
arithmetic’s  base  B  and  the  number  of  base  B  digits  in  the  mantissa  of  a  floating¬ 
point  number  were  not  computed. 

In  this  report,  we  present  an  updated  program  that  computes  the  base  B  from 
the  unit  roundoff  error.  The  new  program  also  computes  the  number  of  base  B 
digits  in  the  mantissa  of  a  floating-point  number  and  the  number  of  bits,  namely  k. 
in  a  base  B  digit.  The  updated  program  calculates  the  page  length  of  virtual  memory 
computer  systems  that  organize  data  in  virtual  memory  as  the  VAX  11/780  does.  In 
order  to  enhance  the  program’s  transportability  between  computer  systems,  ad¬ 
justments  have  been  made  to  the  program  code  and  the  algorithm  for  computing  the 
largest  floating-point  base  B  computer  number. 

Theory  and  Algorithms 


Preliminaries 

If  s  is  a  non-zero,  base  B,  floating-point,  t-digit  computer  number,  wc  assume 
its  representation  is 

s  =  ±Bcf,  -p<e<P 

f  =  ,b,b;  .  .  .  b,,  bj£{0,  1,  2 . B-l },  b,#0.  (1) 

where  B  is  a  positive  integer  power  of  2,  namely  2k,  e  is  in  a  certain  integer  range 
with  both  p  and  P  positive  integers  greater  than  t,  and  the  value  of  f  is  given  by 

f  =  Z  bj/B'  .  (2) 

i  =  I 

In  other  words,  s  is  a  real  number  that  can  be  expressed  as  a  normalized  (b,*0)  base 
B  number  using  at  most  t  digits.  In  particular,  if 

if  i  =  1. 

otherwise,  (3) 
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then  f  assumes  its  smallest  value: 


f  =  B1.  (4) 

If  b,  -  B-l  for  i  =  1, 2,  .  .  .  ,  t,  then  f  assumes  its  largest  value: 

f  =  1  -  B-*.  (5) 

We  assume  that  the  internal  computer  representation  of  each  base  B  digit  b,  is 
given  by  k  two  state  devices  dj,  where  each  device  takes  on  the  value  0  or  1 ;  that  is, 

t>,  =  kH'i1  •  •  •  di"). 

k 

=  Zdf»2H,  d'i(G{0,l}.  (6) 

j  =  i  1  1 

In  particular,  if  f  =  2~r  for  r  =  1 , 2 . k,  then 


r-1  k-r 

b,  =  2k-'  =  (0  .  .  .0  1  0~^~0): 
k 

b;  =  0  =  (6~~?~0):,  i  =  2 . t.  (7) 

Therefore,  the  number  of  bits  in  the  fraction  (mantissa)  f  is  kt.  Thus,  for  a  floating¬ 
point  word  of  fixed  length,  say  N  bits  with  kt  bits  in  the  mantissa,  N-kt-1  bits  are 
reserved  for  the  exponent  e  (stored  as  the  excess  quantity  e  +  a,  a>0)  and  the 
remaining  bit  is  reserved  for  the  sign  of  f.  Furthermore,  since  b,  ?  0,  there  can  be  as 
few  as  kt-k  +  1  significant  (effective)  bits  in  the  fraction  f. 

By  (1)  and  (4)  the  smallest  positive  real  number  in  the  computer  number  set  is 

m  =  B-PB-',  (8) 

whereas  the  largest  positive  real  number  is 

M  =  Bp  (1-B*1)  (9) 

by  (1)  and  (5).  In  a  computation,  the  magnitude  of  results  smaller  than  m  produce 
underflow  (usually  returning  a  base  B,  t-digit  representation  for  zero);  the 
magnitude  of  results  larger  than  VI  produce  overflow. 

The  computer  unit  roundoff  error  u,  w  hich  measures  how  closely  a  real  number 
y  can  be  approximated  by  a  computer  number  s,  is  related  to  the  error  in  s  by  the 
inequality 


iv-si  ^  1  s i u  , 


(10) 
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where 


u 


B'-'  (chops) 
0.5B1*1  (rounds). 


(11) 


The  computer  arithmetic  chops  if  all  of  the  mantissa’s  low  order  digits  b|  beyond  the 
t-th  are  simply  discarded  in  the  normalized  base  B,  infinite  digit  representation  of  y; 
rounds  if  b,  is  incremented  by  1  whenever  bl  +  l>2k~l,  before  discarding  the  low 
order  digits. 


Computation  of  the  Base 

In  the  earlier  program,  we  computed  the  unit  roundoff  error  when  B  ( =  2k)  and 
t  are  unknown.  (The  theory  and  algorithm  for  this  is  in  Appendix  A.)  The  following 
Theorem  and  Corollaries  lead  to  an  algorithm  for  computing  the  base  B  from  the 
unit  roundoff  error. 


Theorem.  In  the  interval  (Be_l ,  Be],  where  -p<e<P,  the  floating-point  numbers  are 
uniformly  spaced  with  spacing  Be'1. 

Proof.  Let  Bef  be  a  floating-point  number  that  is  in  the  interval  [Be-',  Be],  such  that 

f=.b,...bt,  B-'<f<l.  (12) 

To  obtain  the  next  floating-point  number  in  the  interval,  we  add  the  base  B  digit  one 
to  the  last  digit  bt  of  f;  that  is,  the  next  floating-point,  base  B,  t-digit  number  in 
IB'-1,  Be]  is 

B'(f+.b/...b,'),  (13) 

where 

k 

ioTT\ 

k-l 

(6TTTo  d2 

and 

•  b/.  .  .  b/=  Z  6/(2^  =  (2k)_l  =  B  ’.  (14) 

i  =  I 

This  completes  the  proof. 

Corollary  1.  The  floating-point,  base  B,  t-digit  numbers  in  the  interval  [B,_l ,  B1]  are 
uniformly  spaced  with  spacing  B°  =  1. 


if  i  =  1,2,  .  .  .  ,  t— 1, 

if  i  =  t, 


Proof.  Immediate. 
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But  the  value  of  every  floating-point,  base  B,  t-digit  number  in  [BM,  B1)  is  an 
integer.  Therefore,  every  integer  in  [BM,  Bl]  can  be  represented  exactly  as  a 
floating-point,  base  B,  t-digit  number,  since  B'-'  can  be.  In  particular,  the  integers 
2LB'-‘  for  L=0,  1,  2,  .  .  .  ,  k,  have  exact  floating-point,  base  B,  t-digit 
representations. 

Corollary  2.  The  floating-point,  base  B,  t-digit  numbers  in  the  interval  [B<,  B1  r ■]  are 
uniformly  spaced  with  spacing  B. 

Proof.  Since  B'  =  B*1  "-1,  set  e  =  t  +  1  in  the  Theorem.  This  completes  the  proof. 

Thus,  although  B1  can  be  represented  exactly  as  a  floating-point,  base  B,  t-digit 
number  B1  +  1  cannot  be.  Then  how  is  B'-t-  1  approximated  in  the  computer  number 
set? 


Consider  that 


where 


B‘  +  1  =  B1  +  1  B~>  +  B  B  > 

=  Bl  +  I  B1  +  B'"-1  (B-‘B-i) 

=  B1  +  l  (B->  +  B-'->) 

-  BI  +  1  (.  b,  .  .  .  b,  b(  +  1), 

_  (  1  for  i=  1,  t  +  1 

~  \  0  for  i  =  2 . t  . 


(15) 


Therefore,  if  the  computer  arithmetic  rounds,  the  base  B,  t-digit,  floating-point 
approximation  to  B'  + 1  is  given  by 


with 


B'©  1  =  B'  +  > .  b,  .  .  .  bM  b* 


(16) 


if  B  =  2, 

if  B  =  2k  for  k>2  . 


On  the  other  hand,  if  the  computer  arithmetic  chops,  then  for  B  =  2k  (k  =  1 , 2,  .  .  .  ) 
the  base  B,  t-digit,  floating-point  approximation  to  the  sum  B1  +  1  is  given  by  (16) 
with  b*=0. 

Hence,  the  computer  arithmetic  will  approximate  the  real  sum  B'  + 1  by  the 
computer  number  B1,  when  the  computer  arithmetic  rounds  if  B  =  2k  for  k>2;  if 
B  =  2,  then  the  sum  is  approximated  by  Bl  +  B.  On  the  other  hand,  if  the  computer 
arithmetic  chops,  it  approximates  the  sum  B'  +  1  by  B1  for  all  values  of  B  that  are 
positive  integer  powers  of  two. 


4 


Therefore,  given  the  unit  roundoff  error  u  in  (11),  since  its  reciprocal  u~>  is  an 
integer  in  the  interval  (B'-1,  B'J,  and 


2<-u-‘  =  B1 


when 


L  =  {  k-i 

we  have  the  following  algorithm. 

Algorithm  B 

Given  the  unit  roundoff  u,  the  effective  precision,  say  NB1T,  and  whether  the 
computer  arithmetic  rounds  or  chops,  the  algorithm  for  computing  B,  the  binary 
length  of  a  base  B  digit  (namely  k)  and  the  number  of  base  B  digits  in  the  mantissa 
(namely  t)  is: 

B1 .  Set  k  to  zero  (written  symbolically  as  k"- 0). 

s-*-u-' 

s,*-s 

B2.  Do  while  (s:  +  1  *  s,). 

s:-*-2s2 
k-*-k  +  1 
End  Do  while 

B3.  B-s,/s 

B4.  If  (machine  rounds  and  B>2) 

B*-2B 
k—k  +  1 

End  if 

B5.  t*~(NBIT  +  k  -  1  )/k 

Largest  Computer  Number 


(chops) 

(rounds). 


(17) 


In  the  earlier  program  (1  ]  in  order  to  compute  the  largest  computer  number 


TR  6421 


we  begin  by  taking  the  reciprocal  of  the  smallest  positive  computer  number 

m  =  BP-'.  (19) 

If  the  reciprocal  does  not  cause  overflow,  the  program  outputs  an  adjusted  value  of 
the  reciprocal  as  an  approximation  to  M.  Otherwise,  we  find  the  reciprocal  of  the 
product  of  m  by  the  smallest  power  of  two  that  does  not  cause  overflow .  This  gives 


M  =  Bp(0.5).  (20) 

Then  we  output  an  adjusted  value  of  this  result  as  an  approximation  to  M. 

The  problem  with  this  algorithm  is  that  it  requires  invoking  a  routine  (not 
available  on  all  systems)  to  test  whether  an  overflow  occurs.  We  have  corrected  this 
in  the  following  way. 

Assume  that  the  relation  between  the  smallest  and  largest  exponents  of  B  is 
given  by 


p  =  P  +  1.  (21) 

This  is  so  on  the  UNIVAC  1 108.  On  some  machines  the  relation  is  p  =  P;  e.g.,  VAX 
1 1  /780.  Then  w  ith  p  =  P  +  1 


m  =  B-p->  =  B-p-; 


and 


itt1  =  Bp  +  2  =  Bp*-  B 


so  that 


(B3m)‘l(l-BI,)B  =  BP(1-BM) 


(22) 


(23) 


(24) 


approximates  M  without  computer  arithmetic  overflow,  provided  the  order  of 
operation  is  performed  from  left  to  right  as  indicated  in  (24).  For  if  one  were  to 
compute  (B-mV’B  when  p  =  P  1,  the  result  Bp  would  cause  computer  overflow  . 

Similarly,  if  p  =  P,  then 

(B-,m)-|(l-BM)B  =  BP-'U-B1-1)  (25) 


so  that  the  approximation  is  off  by  B. 

Algorithm  M 

Thus,  given  m.  B  and  t,  the  algorithm  for  approximating  M  without  computer 
overflow'  is 


6 


TR6421 


M 1 .  Compute  B-m  and  store  it  in  M. 

M2.  Compute  the  reciprocal  of  \1  and  store  the  result  in  M. 

M3.  Compute  and  store  the  product  Mfl-B1'1)  in  M. 

M4.  Compute  and  store  the  product  MB  in  M.  Output  M. 

Page  Length  Computation 

The  idea  behind  the  algorithm  for  determining  page  length  is  based  on  the  fact 
that  local  data  and  program  code  are  stored  by  the  VAX  1 1  7S0  on  separate  pages 
on  secondary  computer  memory,  with  program  code  following  local  data.* 
Therefore,  if  the  number  of  computer  words,  e.g..  I,  required  by  local  data  occupies 
less  than  a  page  of  secondary  memory,  then  the  first  non-zero  location  following  the 
1th  local  data  word  is  a  program  instruction. 

Algorithm  PL 

Thus,  dimensioning  an  integer  array,  e.g.,  IPGLN,  as  having  length  1,  and 
assuming  the  number  of  local  data  items  is  I,  the  algorithm  for  the  page  length  is 

PL1 .  NMBDT  —  I  +  1 

PL2.  J  —  NMBDT 

PL3.  Do  while  (IPGLN  (J)  =  0) 

J  —  J  +  l 
End  Do  while 

PL4.  J-J-l 

The  final  value  of  J  is  the  page  length. 

Program  Description 

The  updated  single  precision  program  that  incorporates  Algorithms  B,  M  and 
PL  is  given  in  Appendix  B.  Algorithm  u  in  Appendix  A  for  computing  the  unit 
roundoff  error,  which  is  encoded  in  the  earlier  program  [1],  is  implemented  in  the 
new  program  without  any  logical  changes. 

Algorithm  PL  for  computing  page  length  is  implemented  in  the  program  for 
systems  like  the  VAX  11/780  with  page  length  as  great  as  2048  words,  where  the 
local  data  occupies  no  more  than  17  single  precision  words  for  the  single  precision 


'Local  data  refers  to  data  only  uithir,  the  program  module  that  defines  them;  e.g..  data  items  in 
DATA  type  statements  and  intermediate  results  calculated  an  J  stored  oy  the  program  module. 
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version  of  the  program,  and  no  more  than  29  single  precision  words  for  the  double 
precision  version  of  the  program. 

Unlike  the  earlier  program  in  [1],  the  new  program  is  written  in  ANSI  66 
FORTRAN  [2]  (with  the  exception  of  the  PRINT  statement  that  is  AN'S  77  [3])  to 
promote  compilabilitv  on  computers  that  have  different  FORTRAN  compilers. 
Furthermore,  the  program  is  designed  so  that  it  can  be  converted  easily  to  double 
precision  —  simply  remove  the  comment  Hag  C  from  column  1  of  comment  lines  5. 
6,  7,  S6  and  87  in  the  program,  and  convert  lines  8,  84  and  85  to  comment  lines. 


Since  the  program  requires  no  input  data,  it  is  convenient  to  use  in  determining 
the  values  of  the  machine-dependent  parameters  discussed  in  this  report. 

The  program  in  Appendix  B  has  compiled  and  executed  successfully  on  the 
UNI  VAC  1108,  VAX  1 1  -  780,  SEL  32-55  and  ITEL  AS/5.*  The  output  for  the 
double  precision  version  of  the  program  on  the  UNIVAC  1108  and  S  AX  11  780 
appear  in  tables  1  and  2  below,  where  NBIT  is  the  effective  bit  precision  of  a 
floating-point  computer  number;  BASE  is  the  computer  arithmetic's  base  B: 
NBDGT  is  the  number  of  base  B  digits  in  the  fraction  of  a  floating-point  computer 
number;  LNGF1  is  the  number  of  bits  in  a  base  B  digit;  and  RNCHOP  is  the 
arithmetic's  rounding  property  —  0  if  the  arithmetic  chops  results,  1  if  it  rounds. 
The  number  of  bits  in  the  mantissa  of  a  floating-point  computer  number  is  the 
product  of  NBDGT  by  LNGFT 

Table  1.  UNIVAC  1108  Double  Precision 

UNIT  ROUNDOFF  ERROR  =  .173472348-17  NR  J  T  •  60  RNCHOP  =  0 


BASE  =  2.  NBDGT  =  <s0  LNUH 


1 


SMALLEST  F.  P.  NUMBER  = 
APPROXIMATE  LARGEST  F.  P. 
PAGE  LENGTH  -  29 


.273134232213-303 
NUMBER  =  . 898846367431 f 308 


Table  2.  VAX  11/780  Double  Precision 

UNIT  ROUNDOFF  ER9J0  =  0 . 1 3» 7 7 7 ? 7 ? t- 1 6  X'STT  =  56  RNCHOP  =  ! 

BASE  s  2.  NBDOI  =  56  LNGU  *  1 

SMALLEST  F.  P.  VUMPiP  =  0 . 29  3  9  7  35g  770t>D- 36 

A?pF0XT"ATr  LARGEST  F.  ?.  *IU“3ER  =  0 . 8  5  07C  5  9  1  7  1 0  2D*  3  8 

PAGE  LrNGTu  =  176 


■  T no  page  iencr'i  computation  Joe*  rto>  v..»rk  on  : iio  ITF.L  XX  5. 
"T  -Utii/c  ti'ca;  ,:a;a  in  nmial  .l'.e:r,.>r\  a-  the  V  X\  Joe.. 


(he  ITLL.  .went  Joe.  not 
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An  Application  in  Program  Transportability 

Consider  transporting  a  computer  program  from  one  computer  to  another, 
whose  execution  depends  on  the  values  of  machine-dependent  floating-point 
arithmetic  parameters  in  the  program.  In  general,  if  the  parameters  are  not  reset  to 
the  new  machine’s  values,  the  transported  program  will  not  execute  properly.  The 
program  in  Appendix  B,  or  its  double  precision  version,  may  be  used  to  determine 
the  new  parameter  values. 

For  example,  consider  the  segment  of  a  double  precision  program  [4]  in  Ap¬ 
pendix  C  for  calculating  Bessel  functions  Js(x)  and  IN(x)  of  real  argument  and 
integer  order,  which  is  available  on  NUSC’s  UNIVAC  1 108.  The  program  contains 
several  machine-dependent  parameters  that  are  explained  in  the  program’s  com¬ 
mentary  and  are  set  to  the  double  precision  values  of  the  UNIVAC  1108  in  the 
DATA  statement  at  line  75  of  the  program.  Trying  to  execute  this  program  on  the 
VAX  will  fail  as  a  result  of  computer  arithmetic  overflow  if  the  DATA  statement 
values  of  the  parameters  have  not  been  modified  to  reflect  the  VAX’s  double 
precision  arithmetic.  These  parameter  values  can  be  changed  to  reflect  the  VAX's 
values  with  the  aid  of  the  program  in  Appendix  B.  First,  run  the  double  precision 
version  of  the  program  in  Appendix  B  on  the  VAX.  Then  from  its  output  (in  table  2) 
and  the  explanation  of  machine-dependent  constants  in  the  Bessel  function 
subroutine  in  Appendix  C,  one  obtains  the  following  DATA  statement  for  the 
VAX: 


DATA  NSIG.NTEN.LARGEX.EXPARG/ 17. 37, 10000.87. DO 

where  EXPARG  is  the  natural  logarithm  of  the  approximation  to  the  largest 
floating-point  number  in  table  2,  rounded  to  the  nearest  integer. 

Conclusion 

The  program  for  computing  machine-dependent  constants  can  be  used  as  a 
benchmark  for  revealing  features  of  floating-point  number  systems  on  binary 
device  computers,  where  the  base  of  the  number  system  is  a  positive  integer  power 
of  two,  the  mantissa  is  a  normalized  fraction  and  the  computer  .arithmetic  has  at 
least  one  guard  digit.  Furthermore,  the  program  can  be  used  as  an  aid  in  trans¬ 
porting  mathematical  software  containing  machine-dependent  parameters  from  one 
computer  to  another.  In  fact,  it  is  possible  to  facilitate  the  transportation  of  such 
software  between  many  different  computers  by  incorporating  in  the  software 
program  instructions  that  automatically  compute  the  appropriate  machine  values 
for  the  parameters. 
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Appendix  A 


Theory  and  Algorithm  for  I’nit  Roundoff: 
Base  B  and  I  Unknown 


By  Theorem  1,  the  spacing  of  floating-point,  base  B,  t -di-ait  computer  numbers 
in  [B".  B]  is  B1*'.  Therefore,  B0-*-  B1-'  is  in  the  computer  number  set;  furthermore, 
the  real  numbers 


B"  +  vB'-*,  v  =  1,2 . B  '( B- 1 ) 

are  in  [Bu,  B)  and  the  base  B,  computer  number  set;  and.  therefore,  so  are  the 
numbers 


B°  +  2"B1-1,  w  =  0,1,2 . k(t-l)-l, 

where 

2»B‘-»  =  2-l,  L  =  k(t-l)-w, 

are  negative  powers  of  two  that  are  tn  the  computer  number  set,  by  (7). 

However, 

V  =  B"  -r  2-tB'-'  =  BB-'  *  B,-‘(2k'l/B) 

=  B(B->  +  2 k  - 1  /  B ' 1 ') 

=  B(.bj  .  .  .  b,  b. .  ,). 

with 

r  k-1 

(0  ...  0  !);  =  1  if  i  =  1, 

k 

(0  ...  0),  =  0  if  i  =  2 . t. 

k-t 

(1  0  ...  0),  =  :k-i  if  i  =  I  +  1. 

V 

is  a  real  number  that  is  not  in  the  computer  number  set  although  Z^B1  is.  The  real 
number  y  is  approximated  by  a  computer  number  s.  namely 


{B(.b,  .  .  .  b.)  (chops) 
Bt.b,  .  .  .  b,  )  (rounds) 
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m 
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b,  -  (0  .  .  .  0  1 ): 

provided  the  computer  has  a  guard  digit  in  its  arithmetic  hardware  registers  tor 
b, ,  ,  Hence 


(  1  (chopped  arithmetic) 

\  1  +  B1  '  (rounded  arithmetic) . 


Thus,  the  algorithm  tor  computing  the  computer  arithmetic's  relative  accuracv 
u,  effective  precision,  and  rounding  propertv  is 


Algorithm  u 


ul.  Compute  I  <®  t(e  -  2  1 )  for 

L  =  1,2  ...  .  until  either  l*£  =  I  or  (l*£)i-£*l 


u2.  If  (l#£)e£#l,  then  u=£  and  the  machine  rounds;  oiheiwise  u  2i  and  the 
machine  chops. 

u3.  The  minimum  number  of  significant  hits  in  the  mantissa  is  the  Iasi  value  ot  I 
computed  in  ul .  namelv  k< t - 1 )  -  I. 
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Appendix  B 

Computer  Program  for  Calculating  Machine-Dependent  Constants 


■Ml  !f-t 

MMIMtlltlMMItllll 

••  t  Ml  1  Ni>  f  .  t  .  I  ■  It  f  t  f  i If- 1 


i  n . 


hi  ii  fi*i  If  In 


!  liVIJi  Il'N  If  Ml  N  -  !  ■ 

Ill  -I  III  I  «■  t  I  I  L.  !  ,N  .’I  • 

•  '.A,  .f  .  i  .'.  f  .  fi.'.  f 

■  'i  *  .  i  ,’t  f  I  •  ■  IAI  f  •  I . f i f 

!  I  31  R  0  •  HA».  -f  i'll 

■  IM'  nut 

nfi  I  T  /l  f- 1  * 


•  h.’.i  f  .iINf  .  I  UU.IIf/I  T  .UN  f  T,\  ,  Ti  IFF  .f-Ni  Hi Jf  . 

3C 


TUI>|  NNI.fil  0.1. !'• 

rwi.  Nf-fM  i  .ii.' 


Mm.  i  •  OIn*  * 
i .  in 


•  *  *:n  '  ERMINE  UNI’ 


1  < 
1  4 

I  I’ 

I I 
l  ’ 
13 
l? 
30 
s  i 

55 

34 


39 

39 

30 
3  I 
33 

33 

34 

31 

34 

j  i 

36 

35 

40 

41 

43 

43 

44 

43 


iOuni  rtf.  effective  no.  or  his.  (MicnI'Ini.  fRnt.»»» 

UN  11  A  ON  ill.  NF  .  I  INF.  1  ’ 


fin  unt :i  <  ( uni  ta. i  t . uvf  ■ . of 
3  uni r  -  hole  *HNl 1 

UNIT  A  -  UNI!  t  iiMr 

ne  1  t  Nfi  i  r  t  . 

Utfr  ”  UN  I  I A  .N  IT 

IF.  1  UNI  TA.GT  .ONE  '  .A.NI'.  MFF  .f  O.IINf  'F,0  TO  6 
CONTINUE 

IF  .  Ii  IFF  .eo.qne  00  m 
00  TJ  h 

Then  MACHINE  lhOFS 
t  UNIT  -  TuO«.lNIT 

FNCHCF  3EF.0 
'.0  TO  9 

ELSE  MACHINE  F  '  lUNt'S 

8  FNCHCF  =  ONE 

9  CONTINUE 

feint  50.  UNIT.  NUT,  RNCMOF 

*«*r'ETEFMINE  BASI  .  NO.  OF  EASE  M  G  t  T  £  IN  mant;££a.  Mr.  1 7  iNfiHH* 

LNGH  *  0 

SAVE  =■  ONE  UNIT 

SAVE  3  =  SAVE 

I'O  I«HILE<  <SAVE3  ♦  ONE  >  .  NE  .  SAVC3  ) 

10  I F  <  < SAVE 3  +  ONE  .00. SAVE  3  > 00  TO  13 

SAVE  3  *  TU04SA V£  3 
LNGH  -  LNGH  -F  1 
00  ID  10 
CON ' INUE 

13  EASE  *•  SAVE3. SAVE 

IF  (  (  FNCHOF  .EG  .ONE  )  .  ANIi.  (  EASE  .  GT  .  Two  )  MOO  TO  IT 
GO  TO  !4 


4  - 

C  THEN 

13  EASE 

-  TW04EASE 

49 

LNGH 

*  lmc-i  ♦  : 

4  - 

■3  TON’ INUE 

14  .M.CT  = 

'Nfi*  ♦  LNGH  - 

1  '  4  LNGH 

FRINT  53 

.  EASE.  NET'CiT  > 

LNGH 

;  « 

I  »*«CCMFUTE 

SMALLEST  F  .  -  . 

NO  .IM 

5h 

i  ) 


■  j 

C.Z 


:z  IF  1  UNIT .NE . ZERO. AND . UNI TA.NE .UN  I T  >00 
oO  TO  35 

THEN  SMALLEST  r .  F.  NO.  not  REAfHEIi 
3?  UNI ’A  =  UNIT 

UNIT  -  HALF*I)F'IT 
GO  Tn  13 

3.SE  FFINT  SMALLEST  F.  F.  NO . 

35  FEINT  53.  UNITA 

i.ONTINUL 

***AFc'FOXIMATE  LARGEST  F,  F.  ,N0.**« 


past: 


TR  6421 


66 

C 

67 

30  UNIT  *  UNI  rrt*B0SE**3 

6B 

UNIT  -  UNE/UNIT 

69 

UNIT  -  UN  I T  *  < ONE  -  8ASE**(1  -  NBDGT > > 

70 

UNIT  -  UN  I T  4BASE 

71 

C 

CONTINUE 

72 

FEINT  60.  UNIT 

73 

c 

74 

c 

♦  ♦♦DETERMINE  PAGE  LENGTHS 

75 

c 

76 

DO  35  J  -  NMBDT,  7050 

77 

IF  !  I FGLN ( J ) .NE . 0 ) GO  TO  40 

78 

35  CONTINUE 

79 

40  J  -  J  -  I 

80 

IF< J.LT.2049)FRINT  70.  J 

81 

50  FORMA  1 ( 1 H  . 22HUN I T  ROUNDOFF  ERROR  -  . 

EI7.9.5X. 7HNBIT  -  .13 

87 

X  5X.VHRNCH0F'  -  .F2.0  //) 

8  3 

57  FORMAT ! IN  . 7MRASE  =  .F4.0.10H  NBDGT 

*  .  I4.9H  l.NGH  =  .  14  - 

84 

55  FORMAT!  1H  .74HSMALLEST  F.  P  .  NUMBER  - 

,E17.9) 

85 

60  FORMAT  <  1 H  «  35HA!  PROXIMATE  LARGEST  F. 

p.  NUMBFR  =  .  E  1  7.9  i 

86 

c 

55  FORMAT! 1H  .24HSMALLEST  F.  R .  NUMBER  = 

,070. 17) 

87 

c 

60  FORMAT! IH  . 35HAFF ROX I MA IE  IARGEST  F. 

R .  NUMBER  =  . D70 .12) 

88 

•>0  FORMAT!  1H  .  1 4HP  AGE  LENGTH  .14) 

89 

5  TOR 

90 

E  ND 

14 
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Appendix  C 


Double  Precision  Program  for  Calculating  Bessel  Function  JN(\>  and  lN(\) 

M.l|.*5  I  RANC(  1  •  .  f<t  sLh  t 

!  SURKUIJI  I NL  HI'.IKIiX.  NB,  1 7E  .  B.  NCAl  C  ) 

f.  (Hit!  IIJUIINI  CAl  f  ill  AFES  Bl  GSfl  MINI' I  IONS  I  AND  .1  Of  M  AI 

3  I  AKCUMENI  AND  lNIIUEK  UK'Iif  K . 

4  I 

6  L 

*  L  f  <11  i.rlnl  ION  III  I.1  At.  [AM  IS  IN  IHF  (All  ING  HfUlIfNCf 

7  l 

8  l  X  DOUBl  E  I  f  t  •  IS  ION  KEAI  ARGIIMI  NT  Flip  WHICH  1*5  DC  1*5 

«  C  AM  in  RE  LAI  GUI  A  TED.  It  I*S  AM  TO  RE  CAl  CM  All'll. 

10  L  ARSixi  MUST  HE  LEGS  THAN  Ext  AEG  (WHIIH  GM  RFIOWi. 

11  C  rib  IN’IGIR  HIE.  1  +  HIGHFSI  OK  HER  Til  RF  CAl  Llil  AT  F  Ii . 

1.’  C  I-  MUSI  HI  t'OSITiVt  . 

13  C.  I.t  INI  EGER  Mil.  Zf  Ml  IF  ,I*S  ARE  TO  HF  CAl  C.IJ1  AT  f  Ii «  1 

1)  I  IF  I*S  ARE  10  RE  CAL  CULA If  R . 

15  C  H  RUIJRl  E  PRECISION  VFCTOR  OF  LENGTH  NH .  NEC  I'  NOT  HF 

16  C  INITIALIZE  BY  USER.  IF  THE  ROUTINE  lERMINATFS 

l?  C  NORMAL  L  f  (NEALE  *  NH  >  .  IT  RETURNS  JCOR  I  )  -  SIJR  - /E  Ml 

18  r  THROUGH  .H  OR  I > -SUB - NH-M INUS-ONE  OF  X  IN  THIS 

19  L  VECTOR 

CO  r.  NCAl  v.  INTEGER  TYRE  .  NEED  NUT  HE  INITIALIZER  BY  USER. 

Cl  BEFORE  USING  THE  RESULTS.  THE  USER  SHOULU  CHECK  THAT 

CC  C  NCAl  C  NR.  I.E.  ALL  OR IlERS  HAVE  REEN  CALCULATFH  TO 

23  C  THE  HESIREH  ACCURACY.  SEE  ERROR  RETURNS  BELOW. 

24  r 

25  C 

Co  C  EXPLANATION  OF  MACHINE  HEF’ENBENT  CONSTANTS 

C 7  C 

28  C  NS  I G  HECIMAl  SIGNIFICANCE  HESIREH.  SHOULH  BE  SET  TO 

2°  C  IFIX<  Al.0610<2)  *NBI  T  +  1).  WHERE  NB IT  IS  THE  NUMBER  OF 

30  C  BITS  IN  THE  MANTISSA  OF  A  ROUBLE  FRECISION  VARIABLE . 

31  C  SETTING  NSIG  LOWER  WILL  RESULT  IN  RE  CREASED  ACCURACY 

32  C  WHILE  SETTING  NSIG  HIGHER  Will.  INCREASE  CFU  TIME 

33  C  WITHOUT  INCREASING  ACCURACY.  THE  TRUNCATION  ERROR 

34  L  IS  LIMITER  TO  M . 5* 10#* -NSIG  FOR  J*S  OF  OR HE R  LFSS 

35  C  THAN  ARGUMENT.  ANH  TO  A  RELATIVE  ERROR  Or  T  FOR  I* 

36  C  ANH  THE  OTHER  .)#S 

57  C  NTEri  1  ARGES  T  INTEGER  K  SUCH  I  HAT  10**K  IS  MACHINE  - 

38  C  REFRESENTABLE  IN  DOUBLE  FRECISION. 

39  f  LARGE X  UPPER  LIMIT  ON  THE  MAGNITUDE  Or  X.  BEAR  IN  MIND 

40  C  THAT  rr  ABS(X)-N.  THEN  AT  l  EAST  N  ITERATIONS  OF  T HF 

41  C  BACKWARD  RECURSION  WIIL  BE  EXtrUTER. 

42  C  ExFAPG  t  At.GEGT  ROUBLE  Fi-.t.l.  IS  I  ON  ARGUMENT  THAT  THE  1  I BRAF.  Y 

47  L  REXF  ROUTINE  CAN  HANDLE. 

44  C 

45  C 

45  C  ERROR  RETURNS 

4  7  C 

48  C 

4  0  C  LET  G  DENOTE  F5  IT  HER  I  OF  J. 

50  F  IN  CASE  Ilf  AN  ERROk.  N|  111 '.Nf.NB.  AND  NOT  All  i'.*<: 

51  f  ARE  CALCULATED  TO  THF  RE 'IRE  I  ACIIIRACY. 

52  C  IF  NCAl  C.l  T.O.  AN  ARGUMENT  IS  GUI  nr  RANGE.  N  B . .  F  .  0 

53  C  IIP  I7F  IS  NFTTHFR  ('  OR  1  OR  1711  AND  ABS  ■  '  .  i  IF  .  I  YEAR  A  . 

5->  C  1*1  THIS  CASE'.  THE  H  VECTIIR  15  NUT  l.AM  III  A  T  F  R  •  AND  NCAL'. 

■  ■:  15  St.  I  TO  MINO'NB.Ol-l  SO  NCAl  c.ne.nd. 

56  i:  NB.GF.NCALC.GT.O  Will  OCCUR  IF  NR.'.MMACM  AND  ABSip- 

57  C  SUB  -MB -OF  X/G-SUB MAGX  OF -X  I  .  1  T  .  I  O**  ■  N  Tl  N  '  2  .•  .  l.F.  NB 

58  7  IS  MUCH  GREATER  THAN  MAGX.  IN  THIS  CASE  ,  B  ( N  >  IS  cAll'U- 

59  r  LAI2R  TO  THE  RES  If- ED  ACCURACY  FIR  N.IF.NCALC.  BUI  FOR 

AO  L  Hi  Al  f  .1  I  .N.l  F  .NB.  PRECISION  1:.  M’5T.  IF  N.GT.NLAIC  AND 

M  C  ABS'B'NCALCI/BlN) I .FO. 10*#  I  .  THEN  ONLY  THE  FIRST  NSIGK 

C  -.IGNTEICuNT  FIGURES  HE  (it  N  >  MAY  BE  TRIISTFI.  IF  THE  USER 
■•>3  C  WITH  '  TO  CAl  CIILAfE  B'MI  TO  HIGHER  ACCURACY.  HF  SHOHt  R  Ml 

6>  5  AN  AS  YMF  TO  TIC  f  OF.MULA  FOR  HIGH  ORITR. 
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i  <-r<.F.  ifsi  .  ffufA.  n  Hf'K,  u  mcc.f  AM,.‘.niN.sii«r  rows. 

I'l  (Vjr.F'OLti.F  faAVE.fSAVf.1.  ./FW),  If  NTH. HIM  F  »  QUART  *  ONF  »  TWO  .  TEN 
3  .  .Kill I  r  .  Ill  JM 

c  i •  1 1>* 1 1  .KiicHnf-.shAi  t  .ft ii;, cm  in. 

MMLNSIFIN  R'NFi 

I'AT.'i  7ERI1.  If  NIH.OUAk!  .HALF  .UNF  .  I  Wl) .  IF  N/O  .  OHO  ,  .  HiO.  .I'SDO.  .5 
1  .  F'O  ?  I  .  I)  1 . 

I'A  I  A  NS  I G  »N  FFN .LARGE A  *  EXRAkG/  IV.  <07.  lOP’OOf  ‘Mi?/ 
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